clc;clear;
tic;
% PHIP-ALTADENA program for 3 spins of acrylic acid

% Pauli matrices
ai=sqrt(-1);

sx=[0 0.5; 0.5 0];                  % matrix of S_x for spin 1/2
sy=[0 -0.5*ai; 0.5*ai 0];           % matrix of S_y for spin 1/2
sz=[0.5 0; 0 -0.5];                 % matrix of S_z for spin 1/2
od=[1 0; 0 1];                      % 2*2 unity matrix

% evaluating the energy levels of the 6 spin system
nspins=3;
matrixsize=2^3;

% setting {x,y,z} components of spin operators of all spins: 1024*1024
% matrices
I1z=zeros(matrixsize,matrixsize);
I2z=zeros(matrixsize,matrixsize);
I3z=zeros(matrixsize,matrixsize);

I1x=zeros(matrixsize,matrixsize);
I2x=zeros(matrixsize,matrixsize);
I3x=zeros(matrixsize,matrixsize);

I1y=zeros(matrixsize,matrixsize);
I2y=zeros(matrixsize,matrixsize);
I3y=zeros(matrixsize,matrixsize);


% setting matrices of scalar products of all spins of interest
pro12=zeros(matrixsize,matrixsize);
pro13=zeros(matrixsize,matrixsize);
pro23=zeros(matrixsize,matrixsize);


Ix=zeros(matrixsize,matrixsize);
Iy=zeros(matrixsize,matrixsize);
Iz=zeros(matrixsize,matrixsize);


I1z=kron(sz,kron(od,od));
I2z=kron(od,kron(sz,od));
I3z=kron(od,kron(od,sz));

I1x=kron(sx,kron(od,od));
I2x=kron(od,kron(sx,od));
I3x=kron(od,kron(od,sx));

I1y=kron(sy,kron(od,od));
I2y=kron(od,kron(sy,od));
I3y=kron(od,kron(od,sy));

p12=I1x*I2x+I1y*I2y+I1z*I2z;
p32=I3x*I2x+I3y*I2y+I3z*I2z;
p13=I1x*I3x+I1y*I3y+I1z*I3z;

Iz=I1z+I2z+I3z;
Ix=I1x+I2x+I3x;
Iy=I1y+I2y+I3y;


% J-couplings in Hz
J12=10.5;
J32=17.2;
J13=1.8;


tp=100.0;   % preparation time of PHIP: should be long that all coherences are washed out
te=0.0;     % free evolution time of PHIP at the preparation field, should be usually set to zero

% all chemical shifts in ppm in the sigma-scale
ch1=-5.8;
ch2=-6.05;
ch3=-6.25;


% Magnetic fields
% FV: time and number of points:
T2=2.0; % T_2 value, which gives the NMR linewidth 

[BotT] = textread('Bott01mT.DAT');  % B(t) profile of field switching
Nfv=2995;                           % number of points in the profile
        
B=10000*BotT(1,2);                  % starting (preparation) field
Bdet=70400;                         % final (detection) field

H = zeros(64,64);               % Hamiltonian: 64/64 matrix
    rho= zeros(matrixsize,matrixsize);       
    rho0 = zeros(matrixsize,matrixsize);
    rhot = zeros(matrixsize,matrixsize);
    
    rhos=zeros(4,4);
    rhos(2,2)=1/2;
    rhos(3,3)=1/2;
    rhos(2,3)=-1/2;
    rhos(3,2)=-1/2;
    

    H=B*3.0013e2*(ch1*I1z+ch2*I2z+ch3*I3z)/7.04e4;
    H=H+J12*p12+J32*p32+J13*p13;
    
% setting starting density matrix: matrixsize*matrixsize; spins 1 and 2 are ALWAYS
% coming from para-hydrogen

rho0=kron(rhos,od);
rho0=rho0/(2^(nspins-2));



% calculating the transition matrix 


    

% calculating the eigenvenctors
evec=zeros(matrixsize,matrixsize);
evec1=zeros(matrixsize,matrixsize);
energ=zeros(matrixsize,matrixsize);
HH=zeros(matrixsize,matrixsize);
[evec,energ]=eig(H);    % solving the eigen-problem for the Hamiltonian




evec;
evec1=inv(evec);  
        HH=evec1*H*evec;
        
        rho0=evec1*rho0*evec;   % transforming the DM from zeeman basis to the eigen-basis

%check diagonalization!!!        
HH;

%DM evolution at polarization field: washing out the coherences during the
%preparation period

      for i=1:matrixsize
          for j=1:matrixsize
              
              ex=energ(i,i)-energ(j,j);
if ex==0
    rho(i,j)=rho0(i,j);
else
    rho(i,j)=rho0(i,j)*(exp(-2*pi*ai*ex*tp)-1)/(-2*pi*ai*ex*tp)*exp(-2*pi*ai*ex*te);
end              
          end
      end

rhot=rho;

rhot=evec*rhot*evec1;   % now going back to the Zeeman basis; new DM takes into account the evolution at the preparation field


% DM at polarization field is calculated
% now field variation
    Hi = zeros(matrixsize,matrixsize);


for ifv=1:Nfv-1

    Bint=10000*BotT(ifv,2);
    dt=BotT(ifv+1,1)-BotT(ifv,1);
    
    Hi=Bint*3.0013e2*(ch1*I1z+ch2*I2z+ch3*I3z)/7.04e4;
    Hi=Hi+J12*p12+J32*p32+J13*p13;
       
    
   rhot=expm(-2*pi*ai*Hi*dt)*rhot*expm(2*pi*ai*Hi*dt);  % evolution during switching; Hamiltonian is Hi and time period is dt
    
    
    ifv
end

H0=zeros(matrixsize,matrixsize);    % Hamiltonian at the detection field; needed to calculate the NMR spectrum

    H0=Bdet*3.0013e2*(ch1*I1z+ch2*I2z+ch3*I3z)/7.04e4;
    H0=H0+J12*p12+J32*p32+J13*p13;
       
% Applying pulse sequence

fi=pi/2;            % flip angle in NMR
    
    fx=zeros(matrixsize,matrixsize);
    fxx=zeros(matrixsize,matrixsize);
    
    fx=expm(ai*fi*Iy);
    fxx=expm(-ai*fi*Iy);
    
    [e0,v0]=eig(H0);
    e01=inv(e0);
    
    rhot=e01*rhot*e0;   % DM for PHIP-NMR: converting from zeeman basis to eigen basis
    rhoeb1=e01*Iz*e0;   % DM for equilibrium NMR
    fx=e01*fx*e0;
    fxx=e01*fxx*e0;
    magn=e01*(Ix+ai*Iy)*e0;
    
    AAA=zeros(matrixsize,matrixsize);
    BBB=zeros(matrixsize,matrixsize);
    
    % destroying all coeherences in the DM after the spin system arrives to
    % the detection field
                for i=1:matrixsize
                for j=1:matrixsize
               if (abs(i-j)>0)
                   rhot(i,j)=0;
                   rhoeb1(i,j)=0;
               end
                end
                end
    
ZA=fxx*rhot*fx; % PHIP DM after pulse
ZB=fxx*rhoeb1*fx; % DM for equilibrium after pulse
magnt=transpose(magn);

% Detection: 
                
    for i=1:matrixsize
        for j=1:matrixsize
AAA(i,j)=AAA(i,j)+real(ZA(i,j)*magn(j,i));
BBB(i,j)=BBB(i,j)+real(ZB(i,j)*magn(j,i));
                end
    end

 
    
    Nomeg=10000;    % number of pints in the NMR spectrum
    omegF=-6.4;       % staring and final NMR frequency in the spectrum
    omegL=-5.6;
    
   
spe=zeros(Nomeg,3);
% now calculating PHIP-NMR (spe(iw,2) and ordinary NMR (spe(iw,3))
    
    for iw=1:Nomeg
    
            w=omegL-(iw-1)*(omegL-omegF)/(Nomeg-1);
            spe(iw,1)=w;
            w1=Bdet*w*3.0013e2/7.04e4;
for i=1:matrixsize
for j=1:matrixsize
    
    xx=abs(w1-v0(j,j)+v0(i,i));
    k2=1/T2;
    
   spe(iw,2)=spe(iw,2)+AAA(i,j)*k2/(xx*xx+k2*k2);
   
%   real(1.0/(ai*(w1-v0(j,j)+v0(i,i))+1/T2)/sqrt(2*pi));
      spe(iw,3)=spe(iw,3)+BBB(i,j)*k2/(xx*xx+k2*k2);
      
%      /(ai*(w1-v0(j,j)+v0(i,i))+1/T2)/sqrt(2*pi));
end
end

    end

    % plotting the result    
% Plot PHIP NMR
subplot(2,1,1);
plot(spe(:,1),real(spe(:,2)));
title('PHIP-ALTADENA');

% Plot Ordinary NMR
subplot(2,1,2);
plot(spe(:,1),real(spe(:,3)));
title('Thermal equilibrium');
    
% plot(spe(:,1),real(spe(:,2)),spe(:,1),real(spe(:,3)));

% saving the result
save acrylic_acid-tp_100-te_0.dat spe -ASCII
toc;